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Solving  Sparse  Linear  Systems  with  Sparse  Backward  Error 

by 

M  ArioU*.  J.  W.  Demmel**,  I.  S.  Duff 


ABSTRACT 

When  solving  sparse  linear  systems,  it  is  desirable  to  produce  the  solution  of  a  nearby  sparse 
problem  with  the  same  sparsity  structure.  This  kind  of  backward  stability  helps  guarantee,  for 
example,  that  one  has  solved  a  problem  with  the  same  physical  connectivity  as  the  original 
problem.  Theorems  of  Oettli,  Prager  and  Skeel  show  that  one  step  of  iterative  refinement,  even 
with  single  precision  accimiulation  of  residuals,  guarantees  such  a  small  backward  error  if  the 
final  matrix  is  not  too  ill-conditioned  and  the  solution  components  do  not  vary  too  much  in 
magnitude.  We  incorporate  these  results  into  the  stopping  criterion  of  the  iterative  refinement 
step  of  a  direct  sparse  matrix  solver  and  verify  by  numerical  experiments  that  the  algorithm 
firequently  stops  after  one  step  of  iterative  refinement  with  a  componentwise  relative  backward 
error  at  the  level  of  the  machine  precision.  Furthermore,  calculating  this  stopping  criterion  is 
very  inexpensive.  We  also  discuss  a  condition  estimator  corresponding  to  this  new  backward 
error  which  provides  an  error  estimate  for  the  computed  solution.  This  error  estimate  is 
generally  tighter  than  estimates  provided  by  standard  condition  estimators.  We  also  consider  the 
effects  of  using  a  drop  tolerance  during  the  LU  decomposition. 
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1    Introduction 

When  solving  systems  of  n  linear  equations  Ax=b  by  means  of  Gaussian  eliminanon  with  pivoting,  a 
classical  analysis,  (Wilkinson  1961),  shows  that  we  should  expect  to  get  the  exact  solution  x  of  a  slightly 
different  linear  system  ( A+ JA)ii  =  b  +  5b  where  JA  and  5b  are  both  small  with  respect  to  A  and  b.  By 
small  we  mean  small  in  norm,  i.e.  ||5A||  ^AeOAU  and  l|5b|l  ^-fc£||b||  where  ||  .  |1  is  a  matrix  norm,  e  is 
the  machine  precision  (that  is,  the  greatest  positive  number  such  that  fl{\  +  e),  the  floating-point 
representation  of  ( 1  +  £),  equals  1)  and  i  is  the  product  of  the  pivot  growth  factor  and  a  modestly  growing 
fiinction  of  the  dimension  n.  This  classical  view  permits  any  entry  of  5\.  or  Jb  to  be  equally  large,  and  in 
particular  A-i-5  A  may  be  dense  even  if  A  is  quite  sparse.  This  is  unsatisfactory  because  zero  entries  of  A  may 
represent  nonexistent  physical  connections  in  a  system  being  modelled,  and  so  may  be  known  exacdy. 

A  more  satisfying  approach  to  backward  error  than  merely  bounding  ||  5X  \\  and  1|  5h  \\  would  permit  the 
user  to  specify  scaling  factors  eij'>Qandf^  >0  for  each  entry  of  5  A  and  5b,  and  would  compute  the  smallest 
0)>Q  such  that 

\Sa^\<(oeij,\Sb;\<.(of,.  (1) 

By  setting  some  «„•  to  zero,  we  can  insist  that,  if  <b<oo,  the  corresponding  a-  are  known  exactly.  For 
example,  if  «„  =  la;.!  and/^  =  |6;|,  <»  bounds  the  relative  perturbation  in  each  component  of  A  and  b  needed  to 
make  x  an  exact  solution,  and,  in  particular,  5 A  and  5b  have  the  same  sparsity  structures  as  A  and  b.  We 
will  call  this  co  the  componentwise  relative  backward  error.  It  is  important  to  use  this  different  error  estimate 
when  considering  these  restricted  perturbations,  since  Gear  (1975)  has  shown  that  the  conventional  error 
bounds  are  not  appropriate  in  this  case.  It  turns  out  that  the  backward  error  m  is  quite  easy  to  compute,  and  in 
fact  costs  as  little  as  two  matrix-vector  multiplications. 

In  the  following,  if  u  and  v  are  vectors  of  entries  u^  and  V;  and  Q  and  P  are  matrices  of  entries  q^j  and  p  .y , 
Iu|  is  the  vector  of  entries  |  uj,  |Q|  is  the  matrix  of  entries  |  q^\.  u  <  v  means  u^  <  V;  for  all  i,  and  Q  <  P  means 
q^j  <Pij  for  all  i  and;'. 

Theorem  1:  [Oettli  and  Prager  1964]  .  The  smallest  co  satisfying  (1)  is  given  by 

|Ax-b|. 
©=max — :z —  .  (2) 

In  this  expression,  0/0  should  be  interpreted  as  0  and  f/0  (C'^O)  as  infinity.  (o=°o  implies  that  no  co 
satisfying  (1)  exists.  In  particular,  the  smallest  componentwise  relative  perturbation  of  A  and  b  that  makes  x 
an  exact  solution  is 

IAii-b|.- 
6J=max :^ .  (3) 

Thus,  this  theorem  gives  an  a  posteriori  measure  of  the  backward  error  that  is  cheap  to  compute. 

Gaussian  elimination  with  pivoting  does  not  guarantee  that  the  backward  error  co  will  be  small  for  all 
possible  E  and  f.  However,  a  theorem  of  Skeel  (1980)  shows  that  as  long  as  A  is  not  too  LU-conditioned,  and 
as  long  as  the  quantities  (|A|  lx|) ,  in  the  denominator  of  (3)  do  not  vary  too  much  in  magnitude,  then  one  step 
of  iterative  refinement  is  enough  to  guarantee  that  co  will  be  small  for  the  componentwise  relative  backward 
error  in  (3).  This  is  true  even  if  the  residual  r= Ax-b  is  computed  in  the  same  arithmetic  precision  as  used 
for  the  Gaussian  elimination.  The  actual  conditions  under  which  the  foU  jwing  theorem  is  true  are  quite 
complicated,  and  we  refer  for  details  to  Skeel  (1980,  Theorem  5.1) 


Theorem  2:  [Skeel  1980]  Let  £  be  the  machine  precision,  and  let  the  arithmetic  be  such  that  the 
floating-point  result  j7(a 0 6)  of  the  op^ation  a06,(0e  {  +  ,-,x,/})  satisfies /7(aO 6)  =  (aOd)(l+e),  with 
|e|<£.  There  is  a  function  f(.A,b),  typically  behaving  as  0(n),  such  that  when  the  product  of 
K(A)=lllA||A"Mll   and  <T(Ajc)  =  niax(|A|[x|)i/inin(|A||x|)i  is  less  than  (J{A,h)er\  and  there  is  no 

overflow  or  underflow,  the  following  iterative  refinement  algorithm  will  converge  after  one  update  of  x: 
S.olve  Ax=b  using  Gaussian  elimination,  obtaining  solution  x  and  saving  the  LU  factors; 
Compute  the  residual  r= Ax-b  (using  arithmetic  of  machine  precision  s); 
while       <o=max|rj|/(|A||xH-lb|)i>(n+l)£      do 

i 

begin 

Solve  Ad=r  for  d  using  the  saved  LU  factors  of  A; 

Update  x=x-d; 
end; 

This  theorem  may  also  be  extended  to  take  into  account  underflow  and  the  possibility  that,  for  lack  of  a 
guard  digit  in  the  hardware,  we  can  only  assert  that 

fl(a±b)=ail+e^)±b(l+€2). 

where  ICj I  <  £,  (Demmel  1984). 

For  sparse  systems,  it  is  also  possible  to  improve  the  stopping  criterion  of  Theorem  2  by  changing  nto  y, 
the  maximum  number  of  nonzero  entries  in  one  row  of  A. 

Note  that  this  theorem  contradicts  the  usual  advice  that  iterative  refinement  is  not  worth  doing  unless  the 
residual  r=  Ax-b  is  computed  using  arithmetic  of  machine  precision  e^.  Note  also  that  the  theorem  does  not 
say  that  the  refined  solution  will  be  more  accurate,  just  that  it  reflects  the  structure  of  the  original  problem 
more  closely  than  the  unrefined  solution.  If  each  of  the  nonzero  entries  of  the  original  A  is  uncertain  in  its 
least  significant  bit  and  if  w  e,  then  one  could  say  that  one  has  computed  the  solution  as  accurately  as  the 
data  warrants,  since  the  answer  is  exact  for  a  problem  indistinguishable  from  the  problem  one  really  wanted 
to  solve. 

To  use  Theorem  2  as  the  basis  of  a  practical  scheme  for  solving  sparse  linear  systems,  some  modifications 
are  necessary.  In  particular,  when  solving  sparse  linear  systems  where  both  A  and  b  are  sparse  (or  b  has 
components  of  widely  varying  magnitude),  it  often  happens  that  the  quantity  0"(A,x)  in  Theorem  2  is  huge. 
and  convergence  does  not  occur.  Therefore,  we  must  make  another  choice  for  f,  taking  less  account  of  the 
smaller  components  &,-.  This  can  be  done  quite  easily  using  a  modification  of  Theorem  1,  and  is  discussed  in 
Section  2.2. 

There  is  a  new  condition  number  corresponding  to  the  new  definition  of  backward  error  in  (1).  In  the  case 
,.  of  E  =  |A|  andf=lb|,  this  condition  number  is  just  ||  |A"' |  |A|  1| .  This  new  condition  number  is  no  larger  than 
the  traditional  condition  number  ||  A"'  ||  1|  A  ] .  In  fact,  it  may  be  much  smaller  than  ||  A"'  B  ||  A II  if  the  rows 
of  A  are  badly  scaled.  Thus,  combining  the  componentwise  relative  backward  error  with  the  new  condition 
number,  we  obtain  bounds  for  the  real  error  which  are  independent  of  row  scaling.  We  discuss  this  fiuiher  in 
Section  2.1. 

It  has  become  common  to  use  inexpensive  estimators  for  the  usual  condition  number  ||  A"'  ||  ||  A  ||  to 


estiniate  a  bound  for  the  error  in  the  computed  solution  of  Ax=b  (Cline  et  al.  1979,  Higham  1987a, 
Dongarra  et  al.  1979).  In  Section  4,  we  present  an  inexpensive  and  accurate  condition  estimator  for  the  new 
condition  number  ||  |A"'  |  |A|  ||  (and  its  variations).  The  new  condition  estimator  is  based  on  recent  work  by 
Hager  (1984)  and  ffigham  (1987). 

Finally,  we  tested  our  algorithm  and  associated  condition  estimator  in  a  modified  version  of  the  sparse 
linear  system  solver  MA28  (Duff  1977)  from  the  Harwell  Subroutine  Library,  which  uses  the  pivotal 
strategy  of  Markowitz  (1957)  and  a  relative  pivot  test 

Iaj^'|>umax|fl^'| 

on  the  elements  a^  of  the  ;fc-th  pivot  row.  Here  u  (the  threshold  parameter)  is  a  preassigned  fiictor,  usually 
set  to  0.1.  MA28  can  also  drop  entries  of  L  and  U  that  fall  below  a  'drop  tolerance'  in  order  to  further 
decrease  fill-in.  The  L  and  U  factors  are  used  to  solve  Ax  =  b  for  x  by  forward  and  back  substitution  in  the 
usual  way,  followed  by  some  steps  of  iterative  refinement  We  report  on  the  details  of  the  experiments  in 
Section  5.  Our  conclusion  is  that  a  stopping  criterion  like  the  one  in  Theorem  2  (but  suitably  modified  as 
discussed  in  Section  5)  is  a  reliable  and  inexpensive  stopping  criterion  for  iterative  refinement,  often 
stopping  after  one  or  no  update  of  x.  When  drop  tolerances  are  used  and  we  have  convergence,  the  rate  of 
convergence  degrades  slighdy  but  is  still  quite  good-  The  new  condition  estimator  of  Section  4  also  proves 
to  be  inexpensive  to  calculate  and  is  an  accurate  estimate  on  our  test  matrices,  usually  providing  good 
accuracy  for  the  cost  of  a  few  forward  and  back  substitutions  with  the  LU  factors  of  A. 

The  rest  of  this  paper  is  organized  as  follows.  Section  2  discusses  the  componentwise  backward  error 
further  and  also  the  conditioning  of  Ax=  b  with  respect  to  this  backward  error  measure.  Section  3  examines 
how  the  statement  of  Theorem  2  must  change  when  either  the  floating-point  arithmetic  has  no  guard  digit 
(such  as  on  the  CRAY)  or  underflow  occurs.  Section  4  presents  a  condition  estimator  corresponding  to 
componentwise  relative  backward  enor.  Section  5  discusses  the  numerical  experiments.  Section  6  has 
conclusions. 


2    Backward  error  and  conditioning 

2.1    Condition  number 

The  condition  number  of  a  problem  is  the  least  upper  bound  of  the  ratio  of  the  norm  of  perturbation  in  the 
solution  to  the  norm  of  the  perturbation  in  the  input  data,  in  the  limit  as  the  perturbation  in  the  input  data 
goes  to  zero.  To  compute  it,  we  need  a  norm  for  the  perturbation  4  x  in  the  solution  as  well  as  a  norm  for  the 
perturbations  A  A  and  4  b  in  the  input  data.  The  norm  for  the  input  data  will  depend  on  E  and  f  as  described 
above:  |1  (zl  A,zl  b)  ||  ^j  is  defined  as  the  smallest  cysuch  that  \A\.\<(oE  and  |/lb|<cuf.  For  the  norm  of  the 
output,  we  choose  the  usual  sup  norm  |I  x  |L  =  max  lx;|,  in  order  to  cater  for  zero  components  in  x.  With  this 

i 

notation  we  can  write 

<E_f(A,b)  =  Umsup— — -— — —  (4) 

aA-o   ||(4A,4b)|lEr 

4b->0 

where  x+4x=(A+J  A)"'(b+4b).  Following  Skeel  (1979),  this  may  be  easily  evaluated  as 


ll|A-'|ElxMA-MfIL      ■ 

''^•'^^'"^^ m: •  ^^ 

For  example,  if  we  chcx)se  E  =  |A|  and  f  =  |b|  for  the  componentwise  relative  error, 

ll|A-'|lA|txMA-'||b|lL 
'flAUb|(A,b)= jj^^|- (6) 

Sometimes  it  is  convenient  to  have  a  condition  number  which  is  independent  of  the  right-hand  side  b.  Since 

n IA-' I |A|  1x1  i_              , ,  ^, ^^  II IA-' I |A|  tx| IL  ^ 
^^^^^-— ^r^^(A.b)<2 j^^j-— .                                          (7) 

and  II  |A"'  |  |A|  tx|  |I_  / 11  x  IL  ^  D  !A"'  1  |A|  |I„.  we  get  the  simpler  condition  number 

K-^,(A)=lllA-M!A||L>0.5Kr|A,j5,(A,b).  _     (8) 

The  purpose  of  the  condition  number  is,  of  course,  to  provide  error  bounds:  if  A  is  perturbed  by 
|4A|<  6J|A|  and  b  by  |4b|<  a>|bl,  and  if  <a  is  small  enough,  then  x  will  be  perturbed  by  no  more  than  about 
(o  K"|A^|b|(A,b).  More  rigorously.  Skeel  (1979)  shows  that,  for  o)  defined  as  in  (3), 

|5x|L  _aJK-|^|jb)(A,b) 


.<■ 


IxIL       l-a}V|^(A) 


(9) 


Similarly,  if  we  define 

we  have,  for  o)  defined  as  in  (2), 


\  =  in  A  -1 


x:e(A)=I1|A-'1E1L,  (10) 


IJXIL  ^fl)K-Ej(A,b) 


|x|I_      l-(aK-E(A) 


It  is  easy  to  see  thai  the  problem  is  no  more  badly  conditioned  with  respect  to  the  componentwise  relative 
backward  error  measure  than  with  respect  to  the  usual  normed  backward  error  measure.  This  is  because 

x:(A)=  IIA-'  IL  IIA|L>  |11A-'||A||L  =  V|A|(A)  .  (12) 

It  is  possible  for  'r|A,(A)  to  be  much  smaller  than  Ar(A).  For  example,  we  can  make  jc(A)  arbitrarily  large  by 
multiplying  one  of  the  rows  of  A  by  a  large  enough  constanL  However,  k'|^(A)  is  independent  of  the  row 
scaling  of  A. 

2^    Backward  error 

As  stated  in  the  introduction,  it  is  in  practice  necessary  to  modify  the  choice  f  =  |b|  of  the  componentwise 
relative  backward  error.  This  need  arises  because  of  the  factor  a(A,x)  in  Theorem  2;  when  cr(A_x)  is  large, 
convergence  of  the  backward  error  to  in  equation  (3)  to  the  roundoff  level  is  not  guaranteed.  Take,  for 

example,  A  sparse  and  irreducible,  and  x  sparse  such  that  some  6,  =  S  a,yX^  are  zero  because  each  a-Xj  =  0. 

J 
Since  A"'  is  structurally  full  (Duff,  Erisman,  Gear,  and  Reid  1985),  x  will  be  structurally  full  as  well,  so  that 
a  computed  component  i^  can  be  zero  only  through  exact  cancellation.  In  practice,  tiiis  means  that  all 
components  of  the  computed  solution  x  will  be  nonzero,  with  the  entries  which  should  be  zero  containing 
roundoff  error  of  unpredictable  sign.  Therefore  both  /•,  =  (Ax-b)i  and  (|A||x|+|b|).  may  be  small  but  of 
similar  orders  of  magnitude,  so  that  co  stays  large  even  after  some  steps  of  iterative  refinement. 


Ideally,  we  would  like  to  choose  f  to  satisfy  the  following  four  criteria; 

(i)  the  backward  error  co  (in  (2))  usually  converges  to  machine  precision  after  one  step  of  iterative 
refinement, 

(ii)  a)  f  is  "small"  compared  to  b, 

(iif)  the  resulting  error  bound  in  (11)  is  as  small  as  possible,  and 

(iv)  0)  is  row-scaling  independent. 

We  have  experimented  with  two  choices  for  f  which  come  close  to  meeting  these  four  criteria;  this  will  be 
borne  out  by  the  numerical  experiments  in  Section  5.  It  turns  out  we  must  sacrifice  the  sparsity  structure  of 
b  in  order  to  guarantee  a  small  backward  error  bound  o)  (criterion  (i)).  A  trivial  way  to  do  this  is  to  set  E  =  0 
and  f=|r|/e=|Ax-b|/£,  whence  5X=0,  5b  =  r  and  co=e.  Of  course  this  is  unsatisfactory  because  5h  =  r 
may  be  much  larger  in  norm  than  b  if  the  system  is  ill-conditioned,  violating  criterion  (ii).  Our  approach  is  to 
keep  E  =  |A|  and  choose/;  larger  than  |6;|  only  if  it  is  necessary  to  keep  o)  small. 

We  will  choose  f  in  an  a  posteriori  way,  letting  it  depend  on  the  computation  as  follows:  Let  w= |A|  lx|+|b| 
be  the  vector  of  denominators  in  equation  (3).  We  then  choose  a  threshold  r-  for  each  w, ,  so  that  when  w->t- 
we  can  use  the  usual  scaling  factor  fi  =  \b^\.  Otherwise,  when  w-<t^,  we  choose  a  larger  /;. 
Correspondingly,  we  divide  the  equations  of  Ax=b  into  two  categories,  those  where  w^>'C-,  and  those 
where  w^  <  T;  .  We  may  assume  without  loss  of  generality  that  the  leading  m  equations  of  A  x  =  b,  which  we 
denote  by  A^'^x^^^=b^\  belong  to  the  first  category,  and  the  remaining  n-m  equations,  A^'x^^'  =b'^', 
belong  to  the  second.  As  stated  above,  we  will  let  fP^  =  lb^'|  in  the  first  category.  There  are  several 
possibilities  for  r^,  but  in  practice  the  following  one  has  worked  welL  Ti  =  1000ne(II  A,  IL  l|x|L+|feJ), 
where  A.  is  die  ith  row  of  A.  Note  that  T,-  is  about  1000  times  larger  than  the  maximum  possible  roundoff 
error  committed  in  computiing  w-,  and  w,  can  only  be  less  than  r,  if  each  product  a^Xj  is  tiny.  We  performed 
other  runs  to  check  the  sensitivity  of  this  choice  and  found  that  a  change  of  say  a  factor  of  ten  (to  100)  could 
occasionally  change  the  number  of  iterations  and  the  error  estimate  but  usually  not  by  much.  We  note, 
however,  thai  this  can  be  viewed  as  a  local  choice  and  could  be  varied  while  performing  iterative 
refinement,  possibly  increasing  it  in  order  to  decrease  O). 

Given  the  vector  t  of  the  thresholds  T;,  we  can  choose  f^^  in  at  least  two  ways.  The  first  way  that  we 
consider  is  as  follows.  We  let  f^^'  =  |A'^'|e  ||x|L,  where  e  is  thecolumn  vector  of  all  ones.  This  corresponds 
to  the  usual  normwise  backward  error,  and  so  the  components  r,  of  the  residual  are  almost  guaranteed  to  be 
small  compared  to  these /P\  insofar  as  Gaussian  elimination  alone  guarantees  a  small  residual  in  the  norm 
sense.  Since  we  have  not  modified  the  definition  of  E,  we  are  further  guaranteed  a  solution  x  which 
preserves  the  spaisity  structure  of  A. 

There  is  a  difficulty  with  this  choice  of  f,  however  we  are  no  longer  guaranteed  that  ||  Sb  IL  is  small 
compared  to  |1  b  |I_.  This  can  only  happen  when  A  is  very  ill-conditioned,  since  ||  A^'  IL  II  x  IL  /  II  b  IL  is  a 
lower  bound  on  the  condition  number  ||  A".'  |L  II A IL  of  A.  We  have  constructed  artificial  examples  where 
this  happens,  but  not  observed  it  in  practice.  There  is  also  the  possibility  that  large  components  in  f  will 
make  the  condition  number  jf|^|_f  (A,b)  too  large  and  so  make  the  error  estimate  co  x"|^|j(A,b)  too  pessimistic, 
but  note  that  this  condition  number  is  still  bounded  by  2  k-|^|(A).  We  may  avoid  this  possibility  as  follows. 
Given  the  two  backward  errors 
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The  advantage  of  this  formulation  is  that  components  of  f  ^^'  may  be  very  large  compared  to  the  components 
of  b^\  causing  (O^^^  very  small  and  k  to  be  correspondingly  large  but  without  affecting  co^  or  k^^  . 
This  formulation  is  tested  in  the  numerical  experiments  in  Section  5. 

A  second  possible  choice  for  f^'  is  to  use  f ^  =  |I  b  IL  «•  This  choice  of  f^^  assures  us  that  a  small 
backward  error  indeed  means  ll  5b  |L  /  II  b  IL  will  be  small,  but  gives  us  less  assurance  that  the  backward 
error  will  converge  to  machine  precision.  We  have  not  seen  it  fail  in  practice.  As  with  the  other  choice  of  f, 
we  can  bound  the  error  using  two  backward  errors  defined  as  in  (13)  and  the  sum  of  their  products  with  two 
condition  numbers  as  in  (15).  Section  5  also  reports  on  numerical  experience  with  this  backward  error 
measure. 

Both  the  previous  choices  for  f  ^^'  can  violate  one  of  the  criteria  (ii)  or  (iv).  The  choice  f^^  =  |A^^  |  e  ll  x  IL 
guarantees  that  O); ,  i=  1,2,  are  row-scaling  independent  (criterion  (iv)),  while  it  can  violate  criterion  (ii).  The 
choice  f^^'  =  II  b  |I„  e  satisfies  criterion  (ii),  but  the  corresponding  OJj  is  row-scaling  dependent.  Both,  as  we 
shall  see,  satisfy  criteria  (i)  and  (iii). 

We  also  see  that  the  bound  depends  on  the  accuracy  with  which  we  can  compute  the  residual  r  and  the 
backwards  error  <u  in  (2).  How  much  can  roundoff  contaminate  the  computed  (d,  especially  when  r= Ax-b 
is  computed  by  an  arithmetic  with  machine  precision  e?  A  standard  error  analysis  shows  that  the  error  in  the 
computed  r,  5t,  is  bounded  by  ( 7+I)  e(|A|  |xf+|b|),  where  y  is  the  maximum  number  of  nonzero  entries  in  a 
row  of  A.  When  E  =  |A|  and  f=|b|,  this  means  that  the  computed  (o  cannot  differ  from  the  true  (o  by  more 
than  about  ±(  7+l)e  which  will  be  within  the  tolerance  of  our  sparse  modification  of  Skeei's  stopping 
criterion  in  Theorem  2.  Since  the  computed  co  is  almost  certainly  at  least  about  ye,  the  final  error  bound 
(o  x"|^|_|b|(A,b),  can  be  low  by  no  more  than  a  factor  of  2.  The  same  is  true  for  OJ,  ,  /  =  1 , 2. 

At  this  point,  one  might  ask  what  choice  of  E  and  f  minimizes  the  resulting  error  estimate  (1 1).  It  is  easy 
to  see  that  any  choice  of  E  and  f  such  that  Elx|+f  is  a  multiple  of  lr|,  say  E  =  0  and  f  =  |r|,  yields  the  minimum 
product  a;  x-Ej(A,b)=  II  |A-M|r||L/ II  x|L.  Since  the  true  error  is  \\5x\\^l  ||x|i_=  ||  A"'rlL/  i|x|L.wesee 
that  the  bound  is  as  tight  as  ignoring  signs  in  r  allows.  For  this  special  choice  of  E  and  f,  we  should  also  add 
( 7-1-1)  e(|A||x|-i-lb|)  to  |r|  since  roundoff  may  lower  the  computed  value  of  lr|  by  the  same  amount  The 
choice  E  =  0  and  f  =  |r|  -i-  ( y-t-l)  £  (| A|  lx|-i-|b|)  yields  a  new  error  bound  of 
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Thus  we  see  that  the  condition  number  f,AUb((A''')  P'^ys  ^  central  role  independent  of  the  notion  of 
backward  error,  just  because  it  reflects  ±e  possible  roundoff  errors  in  the  computed  residual.  Furthermore, 
after  only  a  few  steps  of  iterative  refinement  Theorem  2  guarantees  that,  to  first  order,  the  bound  (9)  will  be 
alxjut  the  same  as  the  bound  (16).  In  our  experiments  we  have  seen  that,  usually,  the  estimates  of  the  real 
error  given  by  (9)  and  (15)  have  the  same  order  of  accuracy  as  the  estimates  obtained  by  the  bound  (16). 

Note  that  if  we  set  e^-  =  0  A IL  and/;  =  D  b  B.,  the  backward  error  of  x  with  respect  to  E  and  f  is  given  by 
l|Ax-b|L/(||A|L  II X II 1+ II  bIL).  It  is  also  easy  to  see  that 

,,,,     l|A-MLI|A|LI|x||i+llA-'|Lllb|L  „„ 

v,,(A,b)= ^^^ (17) 

which  is  within  a  factor  of  2n  of  |I  A"'  |L  II A  |L-  Thus,  this  choice  of  E  and  f,  which  permits  equally  large 
perturiDations  in  all  enffies  of  A  and  b,  gives  essentially  the  same  backward  error  and  condition  number  as 
the  usual  normed  backward  error. 

We  note,  in  conclusion,  that  Skeel's  original  motivation  (Skeel  1979)  was  to  analyze  the  effects  of  row 
and  column  scaling  of  A  on  the  accuracy  and  the  stability  of  the  LU  factorization.  He  concluded  that  the 
optimal  way  to  scale  depended  on  the  solution:  the  colunms  should  be  scaled  (thus  scaling  the  solution 
components)  so  that  the  components  of  the  scaled  solution  are  all  equal  in  magnitude,  and  the  rows  should 
be  scaled  so  each  component  of  |A|  tx|  (x  is  the  solution)  is  equal  in  magnitude.  This  is  unfortunately  hard  to 
use  in  practice  since  it  requires  much  information  about  the  solution.  Fortunately,  one  step  of  iterative 
refinement  tends  to  overcome  the  effects  of  bad  row  scaling,  as  we  have  seen. 

3    Different  models  of  floating-point  arithmetic 

Theorem  2  assumed  that  arithmetic  was  implemented  rather  cleanly,  i.e.  that  the  floating-point  result 
/7(a06)  of  the  operation  a  06,(0  6  {  +  ,-,x,/})  satisfies 

fl(aOb)={aOb){l+e)  (18) 

with  |ej  <  e,  where  e  is  called  the  machine  precision.  This  model  eliminates  both  the  possibility  of  underflow 
as  well  as  machines  like  the  CRAYs,  where  for  lack  of  a  guard  digit  in  the  hardware  we  can  only  assert  that 

fKa±b)  =  a(l+ei)±b(l+e2)  (19) 

where  \e-\<e.  Thus,  when  a  and  b  are  very  close  and  we  are  subtracting,  this  model  permits  a  large  relative 
error  in  the  computed  difference.  For  example,  on  any  CRAY  or  many  CDC  machines,  the  computed 
difference  of  2'  and  the  next  smaller  floating-point  number  is  wrong  by  a  factor  of  2  (see,  Kahan  1981). 

Despite  this  difficulty,  it  is  possible  to  carry  through  the  pnxjf  of  Theorem  2  using  the  weaker  model  (19) 
instead  of  (18)  and  arrive  at  essentially  the  same  conclusion:  one  step  of  iterative  refinement,  even  without 
computing  the  residual  using  arithmetic  of  machine  precision  £^,  is  enough  to  guarantee  a  small 
componentwise  relative  backward  error  as  long  as  the  matrix  is  not  too  ill-conditioned  and  cr(A,x)  is  not  too 
large.  One  might  expect  problems  in  bounding  the  error  in  the  computed  residual /Z(Ax-b),  since  the  result 
might  be  off  by  a  factor  of  2,  but  in  the  analysis  this  potential  error  is  dominated  by  the  error  in  computing 
Ax,  so  the  proof  goes  through.  Similarly,  the  error  in  updating  x-d  is  swamped  by  larger  errors. 


The  other  exception  to  the  model  in  (18)  is  underflow.  The  extension  of  error  analysis  to  include 
underflow  is  discussed  in  some  detail  by  Demmel  (1984),  and  we  just  summarize  the  results  here.  In  place  of 
(18)  we  use  the  model 

fl(aOb)=(aOb)il+e}^v  (20) 

where  \e\<£as  before,  and  v  represents  the  underflow  error.  Let  A  be  the  underflow  threshold,  that  is  the 
smallest  positive,  normalized  floating-poinl  number.  Then,  on  machines  where  computed  quantines  which 
would  be  smaller  tfian  /I  are  replaced  by  zero,  |v|  is  bounded  by  ^  On  machines  with  IEEE  standard 
floating-point  arithmetic  (see  '^^^^^^  1985,  IEEE  1987),  gradual  underflow  lowers  the  bound  on  |v|  to  eA. 

The  statement  of  Theorem  2  must  be  modified  as  follows  to  account  for  underflow.  For  gradual 
underflow,  we  can  say  the  following:  if  the  inputs  A  and  b  and  the  output  i  are  normalized  (that  is,  exceed  A 
in  magnitude),  and  if  the  residuals  are  computed  by  an  arithmetic  of  machine  precision  either  £  or  £" ,  then 
gradual  underflow  can  only  degrade  performance  to  the  level  of  the  residual  computation  using  the 
arithmetic  of  machine  precision  e.  For  conventional  underflow,  the  norms  of  A,  b  and  x  must  exceed  A/  e  for 
this  statement  to  be  true. 

The  use  of  extended  range  and  precision  in  intermediate  computations  does  not  change  these  conclusions. 
Assuming  r  and  d  are  stored  in  the  same  format  as  A,  b  and  x,  underflows  in  r  and  d  have  the  same  potential 
effects  on  performance  as  they  did  when  they  were  not  computed  in  extended  formaL 

We  have  not  yet  considered  the  effect  of  underflow  on  the  rate  of  convergence  of  the  iteration.  There  are 
matrices  for  which  the  iteration  converges  only  if  underflows  do  not  occur,  but  the  matrices  are  so 
ill-conditioned  as  to  make  the  computed  solution  untrustworthy  anyway.  As  long  as  some  entry  of  A  is  large 
enough  (A  for  gradual  underflow  and  A/  e  for  conventional  underflow)  then  underflows  will  have  an  effect 
on  the  convergence  rate  comparable  to  roundoff. 

4    An  estimator  for  V[^n|(A,b) 

In  order  to  estimate  the  accuracy  of  a  computed  solution  of  Ax=b,  two  ingredients  are  needed:  a  bound 
on  the  backward  error  (however  it  is  measured)  and  a  condition  number  with  respect  to  the  choice  of 
backward  error.  As  discussed  in  Section  2.2,  the  product  of  the  two  previous  quantities  provides  an 
approximate  upper  bound  on  the  relative  error  in  the  computed  solution. 

In  the  case  of  the  conventional  norm  wise  backward  error,  the  condition  number  is  essentially  given  by 
x-(A)=  II  A"'  IL  HAIL.  There  has  been  much  work  on  such  estimators  for  k{A)  in  recent  years  (  for 
example  Cline  et  al.  1979,  see  Higham  1987a  for  a  complete  list  of  references),  and  cheap,  reliable 
estimators  are  available  in  standard  software  packages  such  as  LINPACK  (Dongarra  et  al.  1979).  It  is 
natural  to  seek  an  analogous  estimator  for  K-.^,jb,(A,b). 

From  (5)  we  see  that  the  quantity  we  need  to  estimate  is 

|||A-ME|xMA-MflL=  ll|A-'|(E|x|+f)IL.  (21) 

In  place  of  the  true  solution  x,  we  may  use  its  computed  approximation  x.  In  the  case  of  componentwise 
relative  backward  error,  we  may  also  just  use  the  simpler  condition  number  V|a|(A)  which  requires  us  to 
estimate 

IIIA-'||A||L=ll|A-'||A|e|L  (22) 


where  e  is  the  vector  of  all  ones.  Either  way,  we  need  to  be  able  to  estimate 

IIIA-MgIL  (23) 

where  g  is  a  nonnegative  vector  which  is  easy  to  compute  (in  the  above  examples  it  costs  just  one 
matrix-vector  multiply). 

Let  G = diagCg  1  ,._^  J.  Then  g = G  e  and 

-  « IA-' I  g  II.  =  II IA-'  1 G  e  II-  =  II IA-' GI e  IL  =  II  lA"'  G|  IL  =  II  A"'  G IL-  (24) 

II  A"'G  IL  can  be  estimated  by  the  algorithm  of  Eager  (1984)  and  Higham  (1987),  which  estimates  the 
1-norm  (or  infinity-norm)  of  a  n  x  n  matrix  given  the  ability  to  multiply  a  vector  by  both  the  matrix  and  its 
transpose.  We  can  multiply  any  vector  z  by  the  operator  A~'G  by  multiplying  z  by  the  diagonal  matrix  G, 
and  then  solving  A  y  =  G  z  using  die  LU  factorization  of  A.  Multiplying  by  ( A"'  G)  ^  is  equally  easy. 

Our  estimate  of  condition  numbers  V|^jy(A,b)  includes  a  dependence  on  die  calculated  solution.  We  also 
performed  runs  for  different  solutions  ( for  example,  x^  =  i^,  i=\,...,n  )  and  found  litde  sensitivity.  Note  that 
the  experiments  in  Set  1  in  Section  5  give  us  results  close  to  the  upper  bound  of  twice  k^^. 

5   Numerical  experiments 

We  tested  the  stopping  criteria,  the  backward  errors  (13)  and  the  error  bound  (15)  by  modifying  the  sparse 
linear  system  solver  MA28  in  the  Harwell  Subroutine  Library  (Duff  1977).  As  we  mentioned  in  Section  1, 
MA28  can  drop  entries  of  L  and  U  that  fall  below  a  tolerance  (called  drop  lol  in  our  tables)  in  order  to 
further  decrease  QU-in  {droptol=0  corresponds  to  standard  Gaussian  elimination).  The  resulting  L  and  U 
factors  are  then  used  to  solve  Ax  =  b  for  x  by  forward  and  back  substitution  in  the  usual  way,  followed  by 
some  steps  of  iterative  refinement. 

All  tests  were  done  on  an  IBM  3084.  In  single  precision,  the  machine  precision,  e,  is  16"^  =»  10~*.  In 
double  precision,  it  is  16~"  -2x10"'*. 

AU  our  nins  are  on  a  common  set  of  test  matrices  &om  the  Harwell-Boeing  lest  set  (Duff,  Grimes,  and 
Lewis  1987).  Their  names,  number  of  nonzero  entries  and  condition  numbers  x-(A)  and  k'ia|(A)  are  given  in 
Table  1.  The  name  of  cczh  matrix  includes  its  dimensions,  for  example  GRE115  is  115  by  115.  Two 
matrices  are  identified  as  GRE216.  Both  of  these  have  the  same  structure,  but  they  have  quite  different 
numerical  values.  We  also  ran  our  tests  on  some  odier  matrices  from  the  set  and  obtained  results  broadly 
comparable  with  these  displayed. 

For  each  run,  we  chose  the  value  of  the  solution  x  and  then  we  computed  the  right-hand  side  b  by 
multiplying  the  solution  by  the  test  matrix.  All  matrices  have  also  been  scaled  before  computing  the 
right-hand  side,  thus  obtaining  two  test  problems  for  each  matrix.  The  scaling  is  computed  using  the  Harwell 
routine  MC19,  which  makes  the  nonzeros  of  the  scaled  matrix  near  to  unity  by  minimizing  the  sum  of  the 
squares  of  logarithms  of  the  moduli  of  the  nonzeros  (Curtis  and  Reid  1972).  This  scaling  does  not  guarantee 
that  Ar(A)  and  k'|^|(A)  must  decrease  (see  Table  1)  although  on  many  matrices  the  effect  is  very  beneficial, 
particularly  for  the  classical  condition  number.  This  is  particularly  so  for  the  second  GRE2I6  example, 
where,  before  the  scaling,  the  matrix  was  essentially  singular.  Note  in  general  that  many  of  the  maffices  are 
poorly  conditioned,  particularly  before  scaling. 

In  all  the  runs,  the  standard  normwise  backward  error 
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the  condition  number  k"(A)  and  the  error  bound  t]  KiA)  were  computed  and  compared  to  the  other  backward 
errors,  condition  numbers  and  errOT  bounds. 


Noozoos 

Before  scaling 

After 

1C(A) 

xiA)(A) 

r(A) 

»fiAj(A) 

GRE115 

«i 

a93D+02 

0.86I>+02 

0.69EHO4 

0.13I>tO3 

GRE185 

975 

a38EM36 

0.15D+06 

0J9D+O6 

ai4D+06 

GRE216 

812 

0.7JID+O3 

0.7.7,rHO3 

0.20D+03 

0.18nw)3 

GRE216 

812 

a83I>-15 

0^9D+14 

0.56D+08 

0.85D+07 

GRE343 

1310 

0.47I>+03 

OJTEWB 

030D+O3 

0.26D+03 

GRE512 

1976 

a46I>+03 

0J7D+O3 

0.4OD+O3 

0.36D+O3 

GREllO? 

sm*    . 

0.18EHO9 

0.98D+08 

0.77D+10 

0.24D+09 

WEST67 

794 

O.911>+03 

OJlEWn 

OJOD+03 

0.I3D+O3 

WEST132 

413 

O.llD+13 

0.80D+O7 

0.94I>+O4 

0.21D+04 

WEST156 

■va 

ai2I>+32 

038D+O9 

0.91D+12 

0. 150+09 

WEST167 

506 

0.69D+11 

0J2D+06 

0.46D+O4 

0.12D+O4 

WEs-nsi 

2134 

0J3n*W 

0J8D+O5 

0J8D+O6 

0.53  D+04 

WEST479 

1888 

a49I>-12 

0J7D+O7 

0:271>+O6 

0.20D+O5 

WEST497 

1721 

a38I>-12 

0.13D+07 

0.42EM)7 

0.63  D+04 

WEST6S5 

2808 

0.49EH-12 

0J7D+O7 

0.42I>+O6 

0.36D+Q5 

WEST989 

3518 

0.13D+13 

O.lOD+08 

0^8I>+06 

0.52D+05 

WESnSQS 

5414 

ai4D+13 

O.lOD+08 

0.67D+O8 

0.21D+07 

WECT2021 

7310 

07SD+13 

0.21D+08 

0.86D+06 

O.lOD+06 

Table] 

L  Condition  numbers  before  and  after  scaling. 

We  ran  our  tests  with  different  choices  for  the  vectors  x  and  f  defined  in  Section  22  and  different 
right-hand  sides  b.  According  to  these  different  choices,  we  group  the  experiments  into  3  sets.  We  also 
include  some  runs  using  drop  tolerances  (set  4). 

The  main  data  for  our  numerical  experiments  are  presented  in  Tables  A1-A15  in  the  Appendix.  In  this 
section,  we  display  summaries  of  these  results. 

In  all  cases,  the  stopping  criterion  was 

Stop  if  Q}<  £  or  CO  does  not  decrease  by  at  least  a  factor  of2. 

All  the  runs  used  IBM  double  precision,  except  for  the  experiments  in  single  and  mixed  precision  in  set  1. 
This  stopping  criterion  differs  from  that  used  in  Theorem  2  (a;  <  ( n+  1 )  £).  The  value  in  Theorem  2  can  be 
too  large,  especially  for  very  large  and  sparse  matrices,  and  the  iterative  refinement  could  stop  too  early. 
Generally,  our  stopping  criterion  terminates  the  iterative  refinement  with  a  value  of  o)  less  than  £.  If  the 
convergence  is  slow  (for  example,  using  double  precision,  the  second  GRE216  matrix  in  Table  A7  stops 
after  4  iterations  with  <a=0.4x  10""  =2  e),  our  stopping  criterion  recognizes  this  early.  However,  the  final 
value  of  0)  is  still  of  order  e.  Somewhat  surprisingly  we  find  there  is  no  advantage  in  including  a  factor 
(y+  1)  in  our  stopping  criterion.  Indeed,  its  inclusion  would  often  result  in  no  iteradons,  and  there  are  only 
few  occasions  in  sets  1  to  3  where  die  a)<£  criterion  is  not  met.  Note  that,  in  the  runs  in  sets  2  to  4,  w  is 
replaced  by  co^  +  CO2  (as  in  equations  (13H15)  ).  If  we  used  a  similar  condition  on  77,  in  most  of  the 
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examples  we  did  not  perform  any  steps  of  iterative  refinement  because  the  first  solution  satisfied  the 
stopping  criterion,  but,  before  scaling,  the  estimation  of  the  error  ||  Sx  [\J  ||x  II_  given  by  77  v(A)  was  very 
poor  because  of  the  very  large  value  of  at  (A). 

We  discuss  the  experiments  for  each  of  our  four  sets  of  values  in  turn.  In  all  the  following  tables,  the  row 
corresponding  to  T*Ium.  iter."  gives  the  number  of  steps  performed  by  the  iterative  refinement  algorithm  and 
the  row  corresponding  to  "Num.  cases"  gives  the  number  of  examples  for  which  the  iterative  refinement 
performed  that  number  of  iterations.  By  "Error"  we  denote  the  max-ncrm  of  the  difference  between  the 
computed  solution  and  the  actual  solution  used  to  generate  the  right-hand  side,  divided  by  the  max-norm  of 
the  actual  solution. 

In  the  following,  we  denote  by  co^  and  by  ic^,,  i=  1,2,;=  1,2,3,4,  the  componentwise  backward  enors 
defined  by  (13)  and  the  corresponding  condition  numbers  defined  by  (15).  The  superscript  identifies  the  set 
of  tests. 


Setl: 

For  these  tests  we  chose  T,  =0,  so  that  all  equations  belonged  to  category  1.  Thus  the  backwards  error  was 
given  by  a)f  ^  as  defined  in  (13),  the  condition  number  v^''  and  the  error  bound  by  (O^  k^^^^  as  defined  in 
(15).  Because  all  the  equations  belong  to  category  1,  k^^^  =  x-i^uy  (A,b),  and  coP  =0.  The  right-hand  sides  b 
were  chosen  so  that  the  true  solution  x  had  all  components  equal  to  1.  The  drop  tolerance  was  z  ro.  These 
test  were  nm  in  single  precision,  double  precision,  and  mixed  precision  (all  single  precision,  except  for 
double  precision  computation  of  residuals).  The  Tables  A1-A5  in  the  Appendix  are  relative  to  Set  1. 


.  r(A)  (Before  scaling ) 

^'O   K(A)( After  scaling) 


Log 


,0), 


'"  K';^UAft^r  scaling) 


Log,o(r(A)/ict;^') 


mjn 
-1.9 

-0J8 


■wr 

4.1 

1.4 


19 
&5 


Before  scaling 
min  avr  max 

-0.26  3.6  22 


After  scaling 

min  avT  max 

-0.26  0.91  3J 


Table  2.  Summary  of  results  for  the  condition  numbers  of  set  1. 

In  Table  2,  summarizing  the  results  in  Table  Al,  we  observe  that  the  condition  number  ac^''  is  always  less, 
for  both  scaled  and  unsealed  matrices,  than  twice  the  classical  condition  number  ic(A),  as  must  be  the  case 
from  the  theory.  In  some  examples,  k"^'^  is  much  better  than  k(A)  (for  example,  in  the  WESTl-56  example 
before  scaling  v^'^^  <  3.2  x  10"^  <{A.)).  Moreover,  Table  2  shows  that  the  classical  condition  number 
K-(A),  without  any  form  of  scaling,  is  rather  unreliable  as  a  measure  of  the  ill-conditioning  of  the  system. 
Table  3  (summarizing  the  results  in  Tables  A2  and  A3)  reflects  the  previous  considerations,  so  that  the 
estimation  fuf  «•£"  of  the  error  is  generally  quite  tight,  while  tj  x-(A)  can  be  too  pessimistic  before  scaling. 
Note  that  it  is  possible  for  our  bound  to  be  less  tight  than  that  &t)m  the  classical  theory  but,  when  this 
happens  in  the  experiments,  our  bound  is  only  3  times  greater  than  the  classical  one  in  the  worst  case. 
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Before  scaling 

After  scaling 

Num.  iur. 
Num.  cases 

0 
0 

1 

17 

2   2 

1 

0 
1 

1 
16 

2   2 

1 

min 

■vr 

tpOT 

min 

«VT 

max 

Logi0^o>['h 

-18 
-16 

-16 
-16 

-16 
-16 

-17 
-16 

-16 
-16 

-16 
-16 

^«>o  Error 

a78 

4.7 

22 

a93 

iO 

4.1 

<'<' 

a48 
-0J2 

3.2 

2.5 
20 

0.43 
-0.41 

1.4 
0.53 

3J 
3.0 

^«>''    Error 
r/>f(A) 

- 

Table  3.  Summary  of  results  for  set  1. 

Throughout,  our  estimate  of  condition  numbers  K"iA|jb|(A^'''')  includes  a  dependence  on  the  calculated 
solution.  We  also  performed  runs  for  different  solutions  (  for  example  ,  x^  =  i^,  i=\,...,n  )  and  found  litde 
sensitivity.  Note  that  our  choice  of  x  in  Set  1  gives  us  results  close  to  the  upper  bound  of  twice  k^^.  In 
Tables  A4  and  A5,  we  report  the  results  of  the  algorithm  using  single  and  mixed  precision.  Unfortunately, 
the  test  matrices  are  in  many  cases  so  LU-conditioned  that  the  iterative  refinement  diverged,  that  is  ©{" 
increased  after  some  steps  as  in,  for  example,  GRE1107  and  the  second  GRE216  example  in  Table  A4.  In 
practice,  IBM  single  precision  is  too  poor  to  produce  good  results,  and  the  use  of  mixed  precision  does  not 
help.  Note,  however,  that  our  algorithm  still  terminates  after  only  a  few  steps.  In  every  case,  we  tried 
running  the  iterative  refinement  for  twenty  steps  and  in  no  cases  did  we  get  much  improvement  over  the 
results  shown.  Our  algorithm  for  computing  the  condition  numbers  encounters  numerical  difficulties  pardy 
because  of  the  ill-conditioning  of  these  matrices  and  partly  because  we  use  threshold  pivoting  in  the  LU 
factorization.  We  could  have  used  iterative  refinement  in  this  computation,  but  this  would  be  at  variance 
with  our  desire  for  a  cheap  estimator.  Our  feeling  is  that  single  precision  calculations  are  inappropriate  here. 

Set  2: 

For  diese  tests  we  chose  T;  =  1000n£(IIAi^  |_  I|x|I»+|6.|)  and  f^^=|A^'|e  I|xlL,  where  e  is  the  column 
vector  of  all  ones.  This  leads  to  the  backward  emjrs  <aP  and  tu^'  defined  in  (13)  and  the  condition  numbers 
K-^'  and  K^  and  error  bound  a)P  k'^+o)^'^  k'^2  defined  in  (15).  The  right-hand  sides  were  chosen  so  that 
the  true  solution  x  had  every  fifth  entry  equal  to  1  (x,  =X5=x^  =...=  1)  and  the  rest  zero.  The  drop 
tolerance  was  zero.  These  tests  were  done  in  double  precision  only.  Tables  A6  to  A8  show  the  results  of  runs 
on  set  2.  We  present  a  summary  of  these  results  in  Tables  4  and  5. 

We  also  ran  all  the  test  examples  of  set  2  replacing  zero  with  10"'*  in  x  and  obtained  similar  results.  It  is 
necessary  to  emphasize  that,  in  most  of  the  examples  of  set  2,  the  standard  o)  computed  by  (3)  was  very  large 
(sometimes  of  order  1),  so  that  we  would  get  no  useful  information  if  we  use  a  very  large  value  for  r, .  Notice 
that,  in  ail  our  runs,  cup^  is  very  small  compared  with  (of'^ ,  in  agreement  with  our  comments  after  equation 
(15). 

It  may  appear  that  our  error  estimate  is  sometimes  poor,  but  the  relatively  good  solution  obtained  is  really 
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min 

avT 

max 

,          )C(A)  ( Before  scaJing ) 
*">   K(A)  (  After  icaiiitg) 

-1.9 

4.1 

19 

K^^  (Before  icaling) 

Lap           1 

-0.37 

13 

7J0 

'""  x^HAfter  scaiiAg) 

K^( Before  sealing) 

Log,,.—^ 

'°  n^^/i^ier  scaling) 

-0.43 

U 

6.1 

Before  scaling 

After  scaling 

mm             avT 

max 

min 

avr             max 

Log,a(.<WK^^) 

0.45            43 

23 

0.26 

1.5              3.8 

Log,oMA)/K^) 

0.52            43 

23 

OJO 

1.8              5.2 

Table  4.  Summary  of  results  for  the  condition  numbers  for  set  2. 

fortuitous  as  can  be  seen  by  the  results  in  the  Appendix  using  the  same  matrix  but  with  a  different  right-hand 
side  (the  examples  shown  by  the  GREl  107  results  in  Tables  A3  and  A8  and  by  the  second  GRE216  results 
in  Tables  A2  and  AT). 


Before  scaling 

After  scaling 

Num.  iter. 

0 

1 

2   2 

0 

1 

i   2 

Num.  cases 

1 

13 

4 

2 

12 

4 

min 

avr 

max 

min 

avr 

max 

i^Sio('J) 

-23 

-17 

-16 

-17 

-17 

-16 

^SioC^P) 

-16 

-16 

-15 

-16 

-16 

-15 

^loC"?) 

-32 

-27 

-19 

-31 

-28 

-19 

^^^0  Error 

0.65 

4J 

19 

0.97 

2.2 

4.0 

^S'O              Error 

0.58 

1.7 

4.3 

0.50 

1.6 

2.7 

rjK(A) 

-0.17 

2.8 

16 

-0.23 

0.63 

2.4 

'"<'<'  + -"2™  x:i^ 

Table  5.  Summary  of  results  for  set  2. 

Set  3: 

For  these  tests  we  chose  r,  =  lCXX)n£(I|  A^  ||_  l|x|L+|fe,|)  justas  in  Set  2,  and  f ^^' =  |1  b  IL  e,  where  e  is  the 
column  vector  of  all  ones.  This  leads  to  backward  errors  co^^  and  co^^  defined  in  (13)  and  the  condition 
numbers  ic^^  and  ic^  and  error  bound  co^^  v^Vaj^'  k^^  defined  in  (15).  The  right-hand  sides  were 
chosen  so  that  the  true  solution  x  had  every  fifth  entry  equal  to  1  and  the  rest  zero.  The  drop  tolerance  was 
zero.  These  tests  were  done  in  double  precision  only.  The  Tables  A9-A1 1  are  relative  to  Set  3  of  parameters, 
and  we  summarize  these  in  Tables  6  and  7. 
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- 

min                                avT 

max 

JCtA)  (  BtfoTt  jcaliitif ) 
*'<»  K(.A)(Aft*rscalitig) 

-1.9                               4.1 

19 

K^^(B^orejealiHg) 
^'"'  x^lAfttr  jailing) 

-0J7                              U 

7.0 

K^iB^anseaiing) 

-L9                               4.0 

14 

* 

Before  scaliog 

After  scaling 

min             <vr            max 

mm 

avr             max 

L<7j,oOc(A)/K^) 

0.45             4J              23 

0^ 

IJ              3.S 

Loj,o()r(A)/)C<^) 

0.10           0.97             6.4 

0J8 

0.86             16 

Table  6.  Summary  of  results  for  the  condition  numbers  for  set  3. 


Before  scaling 

After  scaling 

Num.  ittr. 

0 

1 

i   2 

0 

1 

2    2 

Num.  cases 

1 

13 

4 

2 

12 

4 

min 

avr 

max 

min 

avr 

max 

LogiolV) 

-23 

-17 

-16 

-17 

-17 

-16 

i^»,o(<»P') 

-16 

-16 

-15 

-16 

-16 

-15 

^SioC"^) 

-30 

-27 

-17 

-31 

-28 

-19 

^»'0  Error 

0.65 

4J 

19 

0.97 

2.2 

4.0 

^'•o              Error 

0.58 

2.2 

7.9 

0.50 

1.6 

Z7 

,                   rjrCA) 

-0.17 

Z3 

11 

-0.23 

0.63 

2.4 

Table  7.  Summary  of  results  for  the  set  3. 

Comparing  Tables  4  and  6  we  observe  that,  while  ic^^  and  ic^  are  usually  quite  close,  k^  can  be  much 
larger  than  v^^\  (for  example,  see  the  WEST156  example  before  scaling,  where  k^  is  10'*  times  larger 
than  k2^  )  and  the  error  estimation  can  be  pessimistic.  Also  note  that,  comparing  line  7  of  Tables  5  and  7, 
this  choice  of  f  does  not  give  as  good  a  bound  as  our  choice  for  f  in  set  2,  although  the  difference  is  minimal 
after  scaling. 

Set  4: 

For  these  tests  we  used  nonzero  drop  tolerances  {droptol=  10"^,  droptol=  10'').  We  changed  r,  from  its 
earUer  value  to  r.  =  1000rt(e+ci-opwZ)(||  A^  ||_  IIx|L+|6,l)  and  used  f *'' =  I A^"' I e il x IL .  where  e  is  the 
column  vector  of  alJ  ones.  The  entries  of  b  and  x  were  chosen  as  in  Set  3.  Double  precision  was  used.  Tables 
A12-A15  are  the  results  of  runs  using  this  set  of  parameters,  and  the  results  are  summanzed  in  Table  9. 
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We  observed,  contrary  to  Zlatev  (1986),  that  little  gain  in  sparsity  was  obtained  (see  for  example  Table 
A15),  wliile  even  moderate  values  of  drop  tolerance  caused  divergence  of  the  iterative  refinemenL  A  drop 
tolerance  strategy  appears  to  work  well  only  on  very  structured  sparse  matrices  such  as  those  resulting  from 
discretizations  of  partial  differential  equations.  We  confirmed  this  with  a  few  test  runs.  See,  for  example,  the 
results  in  Table  8. 


drop  Id 


10- 


10" 


Fm-in  23619  16085  4697 

NucLiier.  2  14  16 

Eiror  0.32D-14       0.25D-14        0.29D-01 


Table  8.  Fill-in,  numbers  of  iteradons  and  error  for  the  five  point  operator  on  a  30  x  30 
grid,  using   X;  =  1 , i=  l,...ji  and  different  values  of  drop  tol. 


drop 

tol.  -  10- 

-5 

drop  tol.  =  10"' 

Num.  iter. 

0 

1 

2    2 

0                  1 

2;   2 

Num.  cases 

2 

6 

10 

2                  1 

15 

min 

avr 

max 

min             avr 

max 

l^glO^V) 

-18 

-16 

-16 

-18             -15 

-4.6 

^S,o«»{") 

— 

^ 

-17 

-^             — 

-^ 

i^«,o(»r) 

-16 

-16 

-15 

-16             -15 

-2.3 

^*«>  Error 

0.66 

2.Z 

3.7 

0.90             2.1 

4.6 

^^'O              Error 

0.66 

1.6 

2.Z 

0.64             1.6 

3.8 

T7xr(A) 

-0.14 

0.64 

2.5 

-0.95           0.45 

2.9 

''"«,r>ri^;+a,i*)r^' 

Table  9.  Summary  of  results  for  set  4.  The  -«>  entries  correspond  to  values  of  a}{*'  =  0. 

Note  that,  in  this  set,  we  nearly  always  have  o)^*^  =  0.  This  corresponds  to  putting  all  of  the  error  into  b,  that 
is  5A  =  0  and  <5b  =  A  x  -  b,  obtaining  the  situation  which  was  discussed  at  the  beginning  of  Section  12.  In 
this  case,  f  does  not  depend  on  b  explicidy,  but  our  bounds  are  still  good.  Note  again  that  our  stopping 
criterion  terminates  after  only  a  few  iteradons  if  the  iteration  diverges.  We  checked  this  divergence  by 
forcing  more  iterations  and  observed  either  oscillation  or  divergence. 

Finally,  Duff,  Erisman,  and  Reid  (1986,  page  276)  described  an  example  of  Gear  (1975)  where  the  error 
matrix  for  minimizing  the  Frobenius  norm  of  the  error  becomes  arbitrarily  large  if  die  perturbations  are 
constrained  to  the  original  pattern.  On  this  example,  after  one  step  of  iterative  refinement,  using  as  a  starting 
point  die  solution 
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,a=10-•^ 


we  can  guarantee  that  the  error  matrix  E  has  the  same  pattern  as  the  original  matrix.  That  is 

\l    0      0     1/ 
with  a)  <  10"^*  ,  5=  10"*.  It  is  interesting  to  notice  that  x:(A)=  1  +  1/ Jand  K-,A| (A) =4. 

6    Conclusions 

We  have  shown  that,  when  the  iterative  refinement  is  converging,  it  is  possible  and  inexpensive  to 
guarantee  solutions  of  sparse  linear  systems  which  are  exact  solutions  of  a  nearby  system  whose  matrix  has 
the  same  sparsity  structure.  Thus  we  have  answered  the  open  problem  posed  by  Duff,  Erisman  and  Reid 
(1986,  page  276)  concerning  obtaining  bounded  perturbations  while  maintaining  sparsity.  If  the  equations 
arise  from  the  discretization  of  a  partial  differential  equation,  then  a  componentwise  tiny  error  should 
indicate  that  the  solution  obtained  is  that  of  a  neighbouring  partial  differential  equation,  a  conclusion  that 
would  not  be  available  if  classical  error  bounds  were  being  used. 

We  have  extended  work  of  Skeel  (1980)  and  Demmel  (1984)  to  include  the  possibility  of  having  sparse 
right-hand  sides  and  solutions  vectors  and  have  shown  that,  although  we  can  not  always  guarantee  the 
solution  to  a  nearby  problem  whose  right-hand  side  sparsity  is  the  same,  we  can  develop  smtable  bounds  for 
perturbations  in  the  right-hand  side. 

We  discuss  methods  of  inexpensively  and  accurately  calculating  a  condition  number  appropriate  to  this 
tighter  backward  error.  This  condition  number  is  not  bigger  than  that  of  WUkdnson  and  can  indeed  be  much 
smaller,  particulariy  if  the  matrix  is  badly  row-scaled.  For  example,  in  set  1,  the  average  of  the  logarithms  of 
the  ratio  of  the  classical  condition  number  before  and  after  scaling  is  4.1,  while  for  the  Skeel  condition 
number  the  corresponding  value  is  1.4. 

We  have  incorporated  our  backward  error  estimator  in  the  iterative  refinement  step  of  a  direct  sparse 
matrix  solver  and  find  that  we  often  require  zero  or  one  step  of  iterative  refinement  to  guarantee  that  the 
computed  solution  is  the  solution  of  a  nearby  system  with  the  same  sparsity  structure  as  the  original  matrix. 
We  also  observe  that  we  do  not  require  any  extra  precision  in  calculating  residuals,  thus  confirming  remarks 
made  by  Skeel  (1980).  Additionally,  when  combined  with  our  condition  number  estimator,  a  good  estimate 
.  of  the  actual  error  is  obtained.  Furthermore,  when  iterative  refinement  diverges,  our  stopping  criterion 
recognizes  this  early. 

We  observed,  contrary  to  Zlatev  (1986),  that  little  gain  in  sparsity  was  obtained  while  even  moderate 
values  of  drop  tolerance  caused  divergence  of  the  iterative  refinement.  A  drop  tolerance  strategy  appears  to 
work  well  only  on  very  structured  sparse  matrices  such  as  those  resulting  from  discretizations  of  partial 
differential  equations. 

In  this  paper,  we  have  been  using  iterative  refinement  to  improve  the  solution  obtained  using  an  LU 
factorization.  We  have  also  considered  the  case  when  our  LU  factorization  can  be  quite  inaccurate  (set  4).  In 
this  case,  one  could  use  other  techniques  including  SOR  and  CG  and  it  is  a  open  question  as  to  how  far  our 
analysis  could  be  continued  to  cover  these  cases. 
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APPENDIX       Tables  of  results  of  numerical  experiments 

In  all  the  following  tables,  the  column  coiresponding  to  "Num.  iter."  gives  the  number  of  steps  performed 
by  the  iterative  refinement  algorithm.  By  "Error"  we  denote  the  max-norm  of  the  difference  between  the 
computed  solution  and  the  actual  solution  used  to  generate  the  right-hand  side,  divided  by  the  max-noim  of 
the  actual  solution. 


After  scaling 

KA) 

<' 

K-(A) 

GREllS 

0.93D+O2 

0.171>+<B 

a69I>+04 

0.26O+O3 

GRE185 

0J8D+O6 

OJOD+06 

0390+06 

O29O+06 

GRF716 

0.28D+O3 

0.44D+O3 

0200+03 

0J5O+03 

GRE216 

0.831>15 

0J8D+14 

0.56O+08 

0.17I>+08 

GRE343 

0.47r>+03 

0.74CW)3 

O30I>+03 

0.51I>+03 

GRE512 

0.46D+O3 

0.73D+03 

040D+03 

0.72O+03 

GRE1107 

0.18D+O9 

O^D+09 

O77D+10 

0.48O+09 

WEST67 

0.91I>03 

0.15I>+03 

0300+03 

O16I>+03 

WEST!  32 

0.11EM.13 

0.12D+08 

0940+04 

0J3O+04 

WEST156 

0.12D+32 

0J8D+O9 

0910+12 

030I>+O9 

WEST!  67 

0.69CM-11 

0.80D+O6 

046O+04 

0180+04 

WEST381 

0.53D+O7 

0.750+05 

0J38O+06 

O85D+04 

WEST479 

0.49D+12 

0.J7D+O7 

OZ/O+06 

0:250+05 

WEST497 

0J8D+12 

0.200+07 

0420+07 

0120+05 

WEST655 

0.49D+12 

0.57EWJ7 

0420+06 

0410+05 

WEST989 

0.13D+13 

0.16D+08 

0J8O+06 

0.70D+O5 

WEST1503 

0.141>tl3 

0.16D+08 

O67I>+08 

0J5O+07 

WEST2021 

O^EM-13 

0J2D+08 

086O+06 

0.12O+06 

Table  A 1.  Set  1.  Condition  numbers  before  and  after  scaling. 

Num. 

iier.                T) 

1JK(A) 

a.f> 

a>{"  K^l            Error 

GREllS                   I 

0.52D-16 

0.48D-14 

0.59D-16 

0100-13         0.79O-15 

GRE185                   1 

0.12D-15 

0.47D-10 

0160-15 

0.48O-10         0160-12 

GRE216                   1 

0.67D-16 

0.19D-13 

0.67D-I6 

0:290-13         026D-15 

GRE216                   1 

0.73D-16 

0.61D-01 

0110-15 

0.64O-02         0.21O-02 

GRE343                   I 

O.lOD-15 

0.47D-13 

OlOD-15 

0.74O-13         O50D-15 

GRE512                   1 

0.83D-I6 

0J8D-13 

0830-16 

0.61O-13         0.26D-I5 

GREI107                 1 

0.93D-16 

0.17D-07 

0110-15 

0.2?.n-07          0.74D-10 

WEST67                  1 

0.49D-16 

0.45D-13 

0890-t6 

O130-13         O240-14 

WEST132                1 

0.93D-17 

0.98D-05 

0150-15 

O18O-08         O18D-09 

WEST156                1 

0.T/D-I8 

0.9OD+13 

0110-15 

0.42O-O7         0.38D-09 

WEST167                1 

0.80D-I6 

0.55D-O5 

OI20-I5 

0.95D-10         0.48D-11 

WESTiSI                2 

0.45D-16 

0.24D-09 

0160-15 

O12D-10         0:23D-U 

WEST479                1 

0.19D-16 

0.94D-05 

017D-15 

0.96O-09         0.42D-I0 

WEST497                1 

0.77D-16 

079D-04 

Ol  10-15 

0220-09         0.23D-10 

WEST655                1 

0.19D-I6 

0.94D-05 

0.21D-15 

0.12O-08         0.54O-10 

WEST989                1 

0.95D-16 

0.13D-03 

0130-15 

0.21O-08         0.17O-O9 

WEST1505              1 

0.93D-16 

0.13D-O3 

OI6D-15 

0J16D-08         0.17D-09 

WES>-1'2021               1 

0.98D-16 

o:r/u.o3 

OI60-15 

0-52O-08         0.88D-10 

Table  a: 

LSetl.  ii  =  l 

,i=l,...,n ,  double  precisioi 

1  before  scaling. 
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Num.  ilex.                   tj 

77  f  (A) 

<»{" 

«{"  <':i 

Error 

GRE115                    1 

0.64E-16 

0.44E-12 

0.83E-16 

0  7.?F,-13 

0.42E-14 

GRE185                    1 

0.62E-16 

0a4E-10 

0.64E-16 

0.18E-10 

0.54E-13 

GRE216                    I 

0J4E-16 

O.llE-13 

0.79E-16 

0:28E-13 

0.13E-14 

GRE216                    1 

0.89E-16 

OJOE-08 

0.93E-I6 

0.16E-08 

0.17E-09 

GRE343                    1 

0.76E-16 

O^E-13 

0.88E-16 

0.45E-13 

O.lOE-14 

GRE512                    1 

0.76E-16 

OJlE-13 

0.93E-16 

0.66E-13 

0.27E-14 

GREllO?                 1 

0J9E-16 

OJOE-06 

O.lOE-15 

0.48E-07 

0.25E-10 

WEST67                  1 

0J5E-16 

O.llE-13 

0.14E-15 

0.21E-13 

0.89E-15 

WEST132                1 

0.28E-16 

0.26E-12 

0.98E-16 

0J3E-12 

0.73E-14 

WEsnse           c 

)              0J7E-16 

0.52E-04 

0.16E-15 

0.48E-07 

0.98E-08 

wEsner            i 

0.29E-16 

0.13E-12 

O.llE-15 

0.20E-12 

0.44E-14 

WEST381 

I              0.15E-15 

0J8E-10 

0.17E-15 

0.15E-U 

0.56E-I2 

WEST479 

I              0J5E-16 

0.94E-n 

0.22E-15 

0^6E-n 

0.12E-12 

WEST497 

1              0.17E-16 

0.70E-10 

O.llE-15 

0.13E-11 

0.26E-12 

WEST655 

1              0.52E-16 

0.22E-10 

0.19E-15 

0.80E-11 

0.19E-12 

WEST989 

I              0.25E-16 

0.15E-10 

0.12E-15 

0.80E-11 

0.33E-I2 

WEST15Q5 

I              0.50E-I6 

034E-08 

0.17E-15 

0.60E-09 

0.82E-10 

WEST2021 

1              0.50E-16 

0.43E-10 

0.18E-15 

O^E-10 

0.19E-12 

Tables 

\3.Setl.  Xi  =  l, 

f=l^Vi ,  double  precision  after  scaling. 

Num. 

iter.                rj 

»T<f(A) 

<i.P> 

<'<' 

Error 

GRE115 

1             0.15E-06 

O.lOE-02 

0.29E-06 

0.7yH-04 

0.18E-04 

GRE185 

2             0J3E-06 

0.13E+00 

0J3E-06 

0.95E-01 

0.4OE-O2 

GRE216 

1             0J6E-06 

0.73E-O4 

0J9E-06 

0.14E-O3 

0.43E-05 

GRE216 

2             0.59E-06 

0J3E+O2 

0.83E-06 

O.llE+02 

0.43E-01 

GRE343 

1             0J9E-06 

0.11E-Q3 

0.42E-06 

0.17E-03 

0.29E-O5 

GRE512 

1             0.74E-06 

OJOE-03 

0.74E-06 

0.42E-03 

0.15E-04 

GRE1107 

4             0.18E-Q5 

0.13E+04 

O.llE-03 

0.13E+03 

0.86E+00 

wests; 

1             0.15E-06 

0.45E-04 

0.46E-06 

0.19E-Q5 

0.97E-05 

WEST132 

1             0.18E-06 

0.17E-02 

0.47E-06 

0.4tE-04 

0.82E-04 

WEST156 

0             0.22E-07 

0.20E-tO5 

0.54E-06 

0.42E+O1 

0.95E+O0 

WEST167 

1             0.84E-07 

0J8E-Q3 

0.41E-06 

0.19E-04 

0.40E-04 

wEirnsi 

1             0.48E-07 

0.19E-01 

0.51E-06 

O.nE-03 

0.23E-02 

WEST479 

1             0.77F-06 

0.61E-01 

0.95E-O6 

0.62E-03 

0.83E-03 

WEST497 

1             0.12E-06 

0.49E+OO 

0.50E-06 

0.15E-03 

0.17E-02 

WEST6S5 

1             0.74E-O7 

OJlE-01 

0.73E-06 

0.78E-a3 

0.7/H-03 

WEST989 

1             O.llE-06 

0.63E-01 

0.49E-06 

0.89E-03 

0.72E-03 

WEST1505 

1             O.llE-06 

0.73E+OI 

0.70E-06 

0.63E-01 

O.lOE+00 

WES-12021 

1             O.llE-06 

0.93E-01 

0.72E-06 

0  7?F-02 

0^6E-03 

Table 

A4.  Setl.   x,  =  I 

,i=l,...,n,s 

ingle  precis 

on  after  scniing. 
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Num.  iicr. 

V 

r/x-CA) 

<' 

«>r'<' 

Error 

GRE115 

1 

0.20E-06 

0.I4E-02 

0.40E-06 

O.lOE-03 

0.57E-O5 

GRE185 

2 

0.26E-06 

O.lOE+00 

0.58E-06 

0.17E+O0 

0.16E-02 

GRE216 

1 

0J3E-06 

0.66E-04 

0.72E-06 

0.25E-Q3 

OJlE-05 

GRE216 

4 

0.16E-06 

0.89E+O1 

O.llE-05 

0.14E+O2 

0.62E-01 

GRE343 

1 

0J3E-06 

0.97E-04 

0.72E-06 

0.r7H-03 

0.26E-05 

GRE512 

2 

0.25E-06 

O.lOE-03 

0.60E-06 

OJlE-03 

0.72E-Q5 

GREUOT 

4 

0.17E-Q5 

0.12E+O4 

0^E-Q3 

O^E+03 

0.84E+00 

WEST67 

1 

OaOE-06 

0.60E-04 

0.51E-06 

0.21E-Q5 

0.86E-05 

WESn32 

I 

0.15E-06 

0.14E-Q2 

0.75E-06 

0.66E-04 

0.13E-03 

WEST156 

1 

OME-CTJ 

0.98E-h04 

0^9E-06 

0.46E+O1 

0.18E+01 

WEST167 

1 

0.12E-06 

0.53E-03 

0.58E-06 

0.28E-04 

0.16E-04 

WES^nSI 

1 

0.17E-06 

0.67E-01 

0.73E-06 

0.16E-03 

OJlE-03 

WEST479 

1 

0.77E-07 

0.21E-01 

0.63E-06 

0.41E-O3 

0.24E-03 

WEST497 

1 

0.12E-06 

0.51E+00 

0.67E-06 

0.20E-O3 

0.2IE-03 

WEST655 

1 

0.74E-O7 

OJlE-01 

0.82E-06 

0.89E-03 

0.69E-03 

WEST989 

1 

0.94E-(T7 

0^5E-01 

0.88E-06 

0.16E-02 

0.65E-03 

WEST1505 

1 

0.12E-06 

0.80E+01 

0.79E-O6 

0.71E-01 

0.12E+00 

WEST2021 

1 

0.99E-a7 

0.86E-01 

0.80E-06 

0.25E-02 

0.15E-03 

Table  A5.  Set  1.  Xi  =  l, 

i=l,..^ ,  mixed  precision  after  scaling. 

Before  scaling 

After  scaling 

»f(A) 

< 

x-(A) 

'^^. 

.- 

GRBU5 

0.93EM)2 

0J3D+O2 

O^I>+02 

0.69D+04 

0.58CW)2 

0.56D+O2 

GRE185 

0J8D+O6 

OJOD+05 

0.54D+05 

039D+06 

0.461>+fl5 

0.52EM)5 

GRE216 

0^8D+O3 

0.90D+O2 

0.82D+O2 

O^OD+03 

O.llD+03 

O.lOEWn 

GRE216 

0.83EH-15 

0J7D+14 

0.48D+13 

0.56D+O8 

0J5D+07 

0J7D+O7 

GRE343 

0.47I>tO3 

O.I6I>+03 

0.13D+O3 

OJOD+03 

O.lOD+03 

O.llD+03 

GRE512 

0.46D+03 

0.14D+O3 

0.14D+O3 

0.40D+03 

0.14D+O3 

0.14D+03 

GREllO? 

0.18I>+09 

0.4OD+O8 

0J1EM)8 

0.77D+10 

0.91D+08 

0.83D+08 

WEST67 

0.91EM)3 

0^4D+02 

0.78I>O2 

OJOD+03 

0.51D+O2 

0.41D+O2 

WBST132 

O.llD+13 

0.260+07 

0.25D+07 

0.94CM)4 

0.61EM53 

0.83D+O3 

WEST156 

0.12I>-32 

0.12I>+O9 

0.13D+O9 

0.91D+I2 

0.28D+O9 

0.54D+O7 

WEST167 

0.691>i-n 

0.45EM)5 

0J5D+O6 

0.46EW)4 

0.86D+O3 

0.4OD+O3 

WEirn81 

OJ3I>+07 

0.16I>+05 

0.63D+O4 

0J8D+06 

0.23  D+04 

0.13D+04 

WEST479 

0.49D+12 

O.I2EM)6 

072D+07 

0.27I>+06 

0.57EM-04 

0J4IXM 

WEST497 

0J8D+12 

0.75D+O6 

0J3D+O6 

0.42D+07 

0.73D+O3 

0.54D+04 

WEST655 

0.49CH-12 

0.661>+O6 

0.14D+07 

0.42EM)6 

0.12D+05 

0.32D+O4 

WEST989 

0.13E>+13 

0.45rH07 

0.47D+O7 

0.58I>06 

0.21  D+05 

O.llEM-05 

WEST1505 

0.I4D+13 

0.49D+07 

0.53CM-07 

0.67D+08 

0.27D+07 

0.17D+05 

WEl)-r2021 

Oa8I>I3 
Table  A6. 

OJOD+07 
Set  2.  Condi 

0.89D+07 

tion  numbe 

0.S6EM-O6 
rs  before  ani 

0.42D+05 
d  after  seal  in 

O.llD+05 
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Num. 

n 

r7f(A) 

<aP' 

<i,f> 

a,P'  K^* 

Error 

iier. 

a,^  K^ 

GRE115 

0J5D-16 

0.32D-14 

0.84D-16 

0.89D-2S 

0:27D-14 

0.71D-15 

GRE185 

0.94D-16 

0J5D-10 

0.19D-15 

0.24D-25 

0.94D-11 

0.14D-12 

GRE216 

0.I2D-16 

0J4D-14 

0J6D-16 

0.64D-27 

O^OD-14 

0.13D-15 

GRE216 

0.51D-16 

0.42D-01 

0.41D-15 

0.25D-26 

0.15D-01 

0.76D-06 

GRE343 

0.14D-16 

0.65D-14 

0.56D-16 

0.74D-26 

0.90D-14 

O.llD-15 

GRE512 

0.25D-16 

O.llD-13 

0.83D-16 

0.27r>-25 

0.12D-13 

0.19D-15 

GREllO; 

0.42D-16 

0.78D-08 

0.20D-15 

0.58D-24 

0.82D-08 

0.83  D- 10 

WEST67 

0.42D-16 

0J8D-I3 

0.16D-15 

0.27D-30 

0.88D-14 

0.1 2D- 14 

WEST132 

0.24D-16 

0.25D-04 

0.13D-15 

0.80D-28 

0J4D-09 

0.16D-10 

WEST156 

0.12D-22 

0.14D+09 

0.86D-16 

0.I5D-31 

0.100-07 

O.lOD-10 

WEST167 

Oi8D-17 

0.19D-06 

0.20D-15 

0.25D-18 

0.92D-11 

0J7D-12 

WES1381 

0.78D-17 

0.41D-10 

0.15D-15 

0.4OD-29 

0.24D-U 

0.29D-12 

WEST479 

0J3D-19 

0.16D-07 

0J3D-15 

0.14D-28 

039D-10 

0.9  ID- 12 

WEST497 

0.12D-17 

0.44D-06 

0.16D-15 

0.28D-28 

0.12D-09 

OJOD-Il 

WEST655 

0.88D-19 

0.43D-07 

0.26D-15 

0.15D-25 

0.17D-09 

0.29D-11 

WEST989 

0.14D-16 

0.19D-04 

0.14D-15 

0.29D-27 

0.61D-09 

0.26D-10 

WESnSOS 

0.23D-16 

OJlD-04 

0.20D-15 

0.67D-27 

0.99D.09 

0.46D-10 

WES'1-2021 

0.19D-16 

0.52D-04 

0.22D-15 

0J2D-27 

O.llD-08 

0J;4D-10 

Table  A7.  Set 2.  x,.  =  l.i 

:=1,6,..,  elsex-  =  0,  before  scaling. 

Nun* 

n 

»7f(A) 

<«P' 

0)^ 

a,ro  ^m+ 

Error 

iter. 

ajro^ro 

GRE115 

1 

0J2E-17 

0.22E-13 

0.96E-16 

0J6E-27 

0.56E-14 

0.29E-15 

GRE185 

1 

0.64E-16 

0.25E-10 

O.llE-15 

0.41E-24 

0^2E-U 

0J7E-13 

GRE216 

2 

0.60E-16 

0.12E-13 

0.15E-15 

O.lOE-28 

0.16E-13 

0.81E-15 

GRE2I6 

1 

0.12E-15 

0.68E-08 

0.14E-15 

0.94E-25 

0.50E-09 

0.77E-10 

GRE343 

1 

0.60E-16 

0.18E-13 

0.22E-15 

0.48E-26 

0.23E-13 

0.67E-15 

GRE512 

1 

0.86E-16 

0J5E-13 

0.22E-15 

0.25E-25 

OJlE-13 

0.67E-15 

GREllO/ 

3 

O.T/H-16 

0J9E-O6 

0.20E-14 

0.18E-22 

0.18E-06 

O.lOE-08 

WEST67 

1 

0.40E-16 

0.12E-13 

0.16E-15 

0.28E-30 

0.79E-14 

0.13E-14 

WEST132 

1 

0.17E-16 

0.16E-12 

0.17E-15 

0.78E-31 

O.UE-12 

0.54E-I4 

WEST156 

0 

0.6  IE- 17 

0J6E-05 

O.lOE-15 

0.14E-29 

OJOE-07 

0.32E-08 

WEST167 

0 

OaiE-16 

0.94E-13 

0.I8E-15 

0.50E-19 

0.16E-12 

0.24E-14 

WEST381 

1 

0J5E-16 

0.13E-10 

0.12E-15 

0J7E-29 

0.27E-12 

0.86E-13 

WEST479 

2 

0J7E-17 

O.IOE-U 

0.16E-15 

0J3E-30 

0.90E-12 

0.2SE-13 

WEST497 

1 

0.52E-17 

0.22E-10 

O.llE-15 

0.13E-30 

0.81E-13 

0.22E-14 

WEST655 

2 

0.13E-16 

0J5E-11 

0.19E-15 

0.60E-29 

0.22E-n 

0.61E-I4 

WEJIV89 

1 

032E-I6 

0.19E-I0 

oj:oe-i5 

0.63  E-29 

0.43E-11 

0.48E-13 

WEST1505 

1 

0J2E-I6 

0.21  E-08 

0.20E-15 

0.36E-28 

0.54E-09 

0.97E-11 

WEST2021 

1 

0J2E-16 

0.27E-10 

0.20E-15 

0.95E-29 

0.85E-U 

0.18E-I3 

Table  A8.  Set  2.  x,  =  l, 

,1=1,6,...  eis«x,  =  0,  after  scaling. 
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GRE115 

GRE185 

GRE216 

GRE216 

GRE343 

GRE512 

GRE1107 

WESTS/ 

WEST132 

WEST156 

WEST167 

WEST381 

WEST479 

WEST497 

WEST655 

WEST989 

WEST1505 

WEST2021 


«c(A) 

0.93I>02 

0J8D+O6 

0^8I>03 

0.83EM-15 

0.47D+O3 

0.46D+O3 

0.18D+O9 

0.91E>+03 

O.llDi-13 

0.12I>*32 

0.69D+I1 

0^3I>+07 

0.49D+12 

0J8D+12 

0.49D+12 

0.13I>13 

0.14D+13 

0i8D+13 

Table  A9. 


Before  scaling 

'^ 

"j 

0J3I>02 

O^OEHOS 

0.90D+02 

0J7D+14 

0.16DVO3 

0.14D+03 

0.40r>+08 

0J4D+O2 

0.26EMI7 

0.I2D+O9 

0.4SEM)5 

0.16D+O5 

0.12I>+06 

0.75D+O6 

0.66I>+06 

0.45E>+07 

0.49D+07 

OJOD+07 


0J8D+O2 

0.93CM)5 

0.84D+02 

0.18DV15 

O.lOD+03 

0.14I>+O3 

0.42I>08 

0A5D*C1 

0J9D+11 

0.44D+25 

0.68EW)9 

Oi9E>+07 

O^D+12 

O.lOEM-12 

0.18D+12 

0.73D+12 

0.ni>+13 

0.14I>13 


x-(A) 

0.69I>+O4 

0J9CNO6 

QaDX>*m 

0J6EW38 

OJOEMn 

0.40D+O3 

0.T7D+10 

OJOEMB 

0.94D+04 

0.91D+12 

0.46D+O4 

0J8D+O6 

O.270+O6 

0.42D+07 

0.42D+O6 

0^8D+06 

0.67D+08 

0.86D+06 


Aha  scaling 

0^8I>02 

0.4^D+O5 

O.llD+03 

035D+07 

O.lOD+03 

0.14D+O3 

0.91I>+O8 

O^lD+02 

0.61D+O3 

O^SEMW 

0.86D+O3 

O^EMM 

0^70+04 

0.73D+O3 

0.12D+O5 

O^lD+05 

OJZTEMT? 

0.42D+O5 


"i 
0:29D+O4 
0.92D+05 
0.82D+O2 
0.19D+O8 
0.850+02 
0.12I>+O3 
Q2lZh-lO 
0^40+02 
0.27tM)4 
O^D+11 
0.15D+04 
030D+05 
0.28I>+05 
0.85D+06 
0.20EW3S 
O.llD+06 
0.17CM36 
0.12D+O6 


Set  3.  Condition  numbers  before  and  after  scaling. 


GRE115 

GRE185 

GRE216 

GRE216 

GRE343 

GRE512 

GREUOr 

WEST67 

WEST132 

WEST156 

WEST167 

WEST381 

WEST479 

WEST497 

WEST655 

WEST989 

WEST1505 

WEST2021 


Num. 
tier. 

1 
1 
1 
4 
1 
1 
2 
1 
1 
1 
0 
1 
3 
1 
3 
1 
1 
1 


0J5D-16 
0.94D-16 
0.12D-16 
0.51D-I6 
0.14D-I6 
0.25D-16 
0.42D-16 
0.42D-16 
0.24D-16 
0.12D-22 
0.28D-17 
0.78D-17 
OJ3D-19 
0.12D-17 
0.88D-19 
0.14D-I6 
0.23D-16 
0.19D-16 


nx-CA) 

032D-14 

0J5D-10 

0J4D-14 

0.42D-01 

0.65D-14 

O.llD-13 

0.78D-O8 

0J8D-13 

0.25D-O4 

O.HEWW 

0.19D-06 

0.4  ID- 10 

0.16D-07 

0.44D.06 

0.43D-07 

0.19D-04 

OJlD-04 

0.52D-O4 


of 


0.84D-I6 
0.19D-15 
0.56D-16 
0.41D-15 
0^6D-16 
0.83D-16 
0.20D-15 
0.16D-15 
0.13D-15 
0.86D-16 
0.20D-15 
0.15D-15 
0J3D-15 
0.16D-15 
0.26D-I5 
0.14D-15 
0.20D-15 
0.22D-15 


0.89D-28 
0.24D-25 
0.64D-27 
0.25D-26 
0.12D-25 
0J4D-25 
0.58D-24 
0.50D-30 
0.80D-28 
0.17D-27 
0.18D-16 
0.40D-29 
0.14D-28 
0.28D-28 
0.15D-25 
0.29D-27 
0.67D-Z7 
0.32D-27 


0.27D-14 
0.94D-11 
0.50D-I4 
0.15D-01 
0.90D-14 
0.12D-13 
0.82D-08 
0.88D-14 
0J4D-09 
0.75D-03 
0.12D-07 
0.24D-11 
0J9D-10 
0.12D-09 
0.17D-09 
0.61D-09 
0.99D-09 
O.llD-OS 


Enor 

0.7  ID- 15 
0.14D-12 
0.13D-15 
0.76D-06 
O.llD-15 
0.19D-15 
0.83D-IO 
0.12D-14 
0.16D-10 
O.lOD-10 
0J7D-12 
0.29D-12 
0.91D-12 
OJOD-U 
0J19D-n 
0.26D-10 
0.46D-10 
0.24D-10 


Table  A 10.  Set  3.  Xi  =  1 ,  i=l,6,..,  eiyeXi  =  0,  before  scaling. 
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Num. 

n 

riK(A) 

<o['^ 

0,1" 

a,pi  K^^ 

Error 

ilcr. 

<o^  K^ 

GREllS 

1 

0J2E-17 

0.22E-13 

0.96E-16 

0J6E-27 

0.56E-14 

0.29E-15 

GRE185 

1 

0.64E-16 

0.25E-I0 

O.llE-15 

0.41E-24 

0.52E-11 

0J7E-13 

GRE216 

2 

0.60E-16 

0.12E-13 

0.15E-15 

0.12E-28 

0.16E-13 

0.81E-15 

GRE216 

1 

0.12E-15 

0.68E-08 

0.14E-15 

0.94E-25 

0.50E-09 

0.77  E- 10 

GRE343 

1 

0.60E-16 

0.18E-13 

0.22E-15 

0.71E-26 

0.23E-13 

0.67E-15 

GRE512 

1 

0.86E-16 

0J5E-13 

0.22E-15 

a31E-25 

OJlE-13 

0.67E-15 

GREllO; 

3 

0.77H-16 

0.59E-O6 

0.20E-14 

0.18E-22 

0.18E-06 

O.lOE-08 

WEST67 

1 

0.40E-16 

0.12E-13 

0.16E-15 

0.57E-30 

0.79E-14 

0.13E-14 

WESn32 

1 

0.17E-16 

ai6E-12 

O.lTE-15 

0.78E-31 

O.llE-12 

0J4E-14 

WEST156 

0 

0.61E-17 

0.56E-05 

O.lOE-15 

0.14E-29 

OJOE-07 

0J2E-08 

WEST167 

0 

0.21E-16 

0.94E-13 

0.18E-15 

O^OE-19 

0.16E-12 

0.24E-14 

WESTJ81 

1 

0J5E-16 

0.13E-10 

0.12E-15 

0.57E-29 

0.27E-12 

0.86E-13 

WEST479 

2 

037E-17 

O.lOE-n 

0.16E-15 

033E-30 

0.90E-12 

0.28E-13 

WEST497 

1 

0.52E-17 

O^E-10 

O.llE-15 

0.13E-30 

0.81E-13 

0.22E-14 

WEST655 

2 

0.13E-16 

0.55E-11 

0.19E-15 

0.60E-29 

0.22E-11 

0.61E-14 

WEST989 

1 

032E-16 

0.19E-10 

0.20E-15 

0.63E-29 

0.43E-11 

0.48E-13 

wEs-nsos 

1 

0J2E-16 

0.21E-08 

OaOE-15 

0J6E-28 

0.54E-O9 

0.97E-1I 

WEST202I 

1 

0J2E-16 

0.27E-10 

0.20E-15 

0.95E-29 

0.85E-11 

0.18E-13 

Table  All. ; 

Set  3.    X.  =  l, 

,  i=1.6,... 

elsex-=0,  after  scaling. 

droptol  =  10" 


droptol  =10" 


«f(A) 

<' 

< 

<' 

GREllS 

0.69E-K}4 

O.OOE+00 

0.12E+O3 

O.OOE+OO 

0.12E+03 

GRE185 

0J9E+06 

O.OOE+00 

0.17E+06 

O.OOE+00 

0.14E+06 

GRF716 

0.20E+O3 

O.OOE+OO 

0.21E+O3 

O.OOE+00 

0J21E+O3 

GRF,716 

0.84E+O8 

O.OOE+OO 

0.15E+O8 

O.OOE+OO 

O.lOE+07 

GRE343 

OJOE+03 

O.OOE+00 

OJlE+03 

O.OOE+OO 

0.26E+03 

GRE512 

0.40E+03 

O.OOE+00 

0.43E+O3 

O.OOE+00 

0J7E+03 

GRE1107 

0.63E+10 

O.OOE+OO 

0.23E+O9 

O.OOE+OO 

0J5E+O7 

WESTOT 

OJOE+03 

0.29E+01 

0.16E+03 

O.OOE+OO 

0.14E+03 

WEST132 

0.94E+O4 

O.OOE+OO 

0.24E+04 

O.OOE+00 

0.77F,+O4 

WEST156 

0.91E+12 

O.OOE+OO 

0.29E+O9 

O.OOE+OO 

0.16E+O6 

WEST167 

0.46E+O4 

O.OOE+OO 

0.16E+04 

O.OOE+00 

0.13E+04 

WES-rJ81 

0J8E+06 

O.OOE+00 

0.d5E+04 

O.OOE+00 

0.54E+04 

WEST479 

0.27E■^06 

O.OOE+OO 

0.23E+05 

O.OOE+00 

0.20E+05 

WEST497 

0.42E+fl7 

O.OOE+OO 

0.65E+04 

O.OOE+OO 

0.63E+O4 

WEST6S5 

0.42E+O6 

O.OOE+OO 

0.43E+O5 

O.OOE+00 

0.37E+05 

WEST989 

038E-^06 

O.OOE+OO 

0.63E+O5 

O.OOE+OO 

0.53E+05 

WEST1505 

0.67E+08 

O.OOE+00 

0J5E+07 

O.OOE+00 

0.21E+O7 

WEi>l-i021 

0.86E+06 

O.OOE+OO 

0.I2E+06 

O.OOE+00 

O.lOE+06 

Table    A12.    Set 

4.    Condition    numbers    after 

scaling    for 

droptol. =  10 

drop 

ro/.=  10"\ 

and 
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Num. 
iler. 

7 

r\K(X) 

«<" 

««(*' 

Error 

GRE115 

2 

0.99E-18 

0.68E-14 

0.OOE+0O 

0.69E-16 

0.85E-14 

0.15E-14 

GRE185 

3 

0.55E-16 

a22E-10 

O.OOE+00 

OJOE-16 

0.83E-11 

0.80E-13 

GRE216 

1 

0.90E-16 

0.18E-13 

O.OOE+00 

0.83E-16 

0.18E-13 

0.88E-15 

GRE2I6 

29 

O.lOE-15 

0.84E-08 

O.OOE+00 

OJOE-15 

0.T7E-08 

0.63E-I0 

GRE343 

1 

0.90E-16 

0.27E-13 

aOOE-eOO 

0.83E-16 

0.26E-13 

0.81E-15 

GRE512 

1 

0.86E-16 

0J5E-13 

O.OOE+00 

O.llE-15 

0.48E-13 

0.68E-15 

GRE1107 

15 

0.62E-16 

0J9E-06 

O.OOE+00 

0.19E-15 

0.44E-07 

a27E-09 

WEST67 

1 

OJOE-16 

0.15E-13 

0.13E-I6 

0.61E-16 

O.lOE-13 

0.85E-15 

WEST132 

2 

0J6E-16 

0J3E-12 

O.OOE+00 

0.67H-16 

0.16E-12 

0.46E-14 

WEST156 

0 

0.61E-17 

0J6E-05 

O.OOE+00 

0.54E-16 

0.16E-07 

0J2E-08 

WEST167 

0 

0.21E-16 

0.94E-13 

O.OOE+00 

0.67E-16 

O.llE-12 

0.24E-14 

WEST381 

2 

0.23E-16 

0.89E-11 

O.OOE+00 

0J4E-16 

036E-12 

0.78E-13 

WEST479 

3 

0.26E-16 

0.71E-11 

O.OOE+00 

0.57E-16 

0.13E-11 

0J5E-13 

WEST497 

1 

0J8E-17 

0.25E-10 

O.OOE+00 

0J5E-16 

0J6E-12 

0.47E-14 

WEST655 

2 

0J5E-16 

0.23E-10 

O.OOE+OO 

0.91E-16 

0J9E-11 

0.22E-I3 

WEST989 

1 

0.13E-15 

0.75E-10 

O.OOE+00 

0.19E-15 

0.12E-10 

0.18E-13 

WEST1505 

2 

0.64E-16 

0.43E-08 

O.OOE+OO 

O.lOE-15 

0J5E-09 

O.lOE-10 

WEST2021 

2 

0.95E-16 

0.82E-10 

O.OOE+00 

0.13E-15 

0.16E-10 

0.59E-I3 

Table  A13. 

Set  4.    Xi  = 

1.1=1,6,... 

elseXi  =  0,  after  scaling  and  drop  tol 

.=  lo-^ 

Nirni, 
iler. 

n 

T/fCA) 

a,r 

«>|*» 

Error 

GREllS 

4 

0J5E-17 

0.24E-13 

O.OOE+00 

0.48E-16 

0.59E-14 

0.80E-15 

GRE185 

15 

0.46E-16 

0.15E-10 

O.OOE+00 

0.61E-16 

0.87E-n 

0.14E-12 

GRE216 

1 

0.65E-16 

0.13E-13 

O.OOE+00 

0.74E-16 

0.16E-13 

O.llE-14 

GRE216 

3 

0.26E-04 

0.15E+O3 

O.OOE+00 

O.llE-02 

0.12E+04 

0.22E+01 

GRE343 

3 

0.66E-16 

0.20E-13 

O.OOE+OO 

0.87E-16 

0.23E-13 

0.72E-I5 

GRE512 

4 

0.63E-16 

0.26E-13 

O.OOE+OO 

0.89E-16 

0J2E-13 

0.79E-15 

GRE1107 

3 

0.64E-05 

O.lOE+04 

O.OOE+OO 

0.16E-02 

0.90E+04 

0.13E+01 

WEST67 

2 

037E-16 

O.llE-13 

O.OOE+OO 

0.45E-16 

0.61E-14 

0.14E-14 

WEST132 

3 

0::5E-16 

0.23E-12 

O.OOE+OO 

0J2E-I6 

O.llE-12 

OaiE-14 

WEST156 

0 

0.59E-18 

0.73E-08 

O.OOE+00 

0J4E-16 

0.87E-U 

0.18E-12 

WEST167 

0 

0.21  E- 16 

0.94E-13 

O.OOE+OO 

0.67E-16 

0.84E-13 

0.24E-14 

WEST381 

4 

0.17E-16 

0.67E-11 

O.OOE+OO 

0.53E-16 

0.29E-12 

0J3E-13 

WEST479 

7 

034E-17 

0.91E-12 

O.OOE+OO 

0.55E-16 

O.llE-11 

0.51E-13 

WEST497 

4 

OJOE-17 

0.13E-10 

O.OOE+OO 

0.58E-16 

0J6E-12 

0J6E-14 

WEST655 

5 

0:28E-16 

0.12E-10 

O.OOE+OO 

0.65E-16 

0.24E-11 

0.55E-13 

WEST989 

5 

0J2E-16 

0.19E-I0 

O.OOE+OO 

0.64E-16 

0.34E-11 

O.llE-12 

WEST1505 

10 

0J2E-16 

0.20E-08 

O.OOE+OO 

0.90E-16 

0.19E-09 

0.23E-10 

WEST2021 

5 

0J2E-16 

o:r/E-io 

O.OOE+OO 

0.94E-16 

0.98E-11 

0.72E-13 

Table  A 14. 

Set  4.   x.= 

=  l.i=l,6,.. 

,  ciyejCj=0,  after  scaling 

and  drop  tol.  =  \Q~^ . 
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Nonzeros  Fill-in 

droptol=Q.O    drop lol=lO~' drop u>l=  10'^ 


GRE115 

421 

647 

651 

605 

GRE185 

975 

3173 

3028 

2929 

GRE216 

812 

2544 

2263 

2262 

GRE216 

812 

2767 

2580 

2180 

GRE343 

1310 

5334 

4891 

4890 

GRE512 

1976 

11535 

11020 

11007 

GREU07 

5664 

47603 

45255 

41181 

WEST67 

294 

267 

202 

204 

west:  32 

413 

89 

87 

83 

WEST156 

362 

27 

20 

15 

WESTlff/ 

506 

96 

96 

92 

WES1381 

2134 

2057 

1867 

1711 

WEST479 

1888 

1121 

982 

790 

WEST497 

1721 

279 

263 

252 

WEST655 

2808 

2092 

1791 

1709 

WEST989 

3518 

1156 

1139 

1135 

WEST1505 

5414 

2032 

1934 

1821 

WEST2021 

7310 

2539 

2466 

2410 

Table  A15.  Set  4.  Number  of  nonzero  entries  in  the  original  matrices  and  fill-in  for 
drop  tol.  =  0.0,  drop  lot.  =  10"^  and  drop  tol.  =  10"'  after  scaling. 
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